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A mesh objective crack band model was implemented within the generalized method 
of cells micromechanics theory. This model was linked to a macroscale finite element 
model to predict post-peak strain softening in composite materials. Although a mesh 
objective theory was implemented at the microscale, it does not preclude pathological 
mesh dependence at the macroscale. To ensure mesh objectivity at both scales, the energy 
density and the energy release rate must be preserved identically across the two scales. 
This requires a consistent characteristic length or localization limiter. The effects of scaling 
(or not scaling) the dimensions of the microscale repeating unit cell (RUC), according to 
the macroscale element size, in a multiscale analysis was investigated using two examples. 
Additionally, the ramifications of the macroscale element shape, compared to the RUC, 
was studied. 


I. Introduction 

An advantage composites possess over most monolithics as a structural engineering material is that, 
in addition to the physical properties of the constituents, the subscale geometry (architecture) of the con- 
stituents contribute to the apparent response of the composite, yielding a material which exhibits behavior 
that surpasses the sum of its constituents. Often, these details are smeared into a homogenized model for 
convenience, and ad-hoc assumptions are used to include the effects of the composite micro-architecture, in 
computational analysis methods. These assumptions work well to predict the elastic behavior of composite 
materials. However, to incorporate the non-linear effects of damage and failure, increasingly complicated 
continuum damage theories and failure criteria must be developed, such as those included in numerous review 
papers and books. 1_ ' Many, but not all, of these theories incorporate various non-physical parameters that 
must be calibrated in order to capture the appropriate failure modes in the composite. 

The deficiencies of homogenized models become more apparent when strain localization leading to soft- 
ening damage occurs in the material. Once localization occurs, the characteristic length of the material 
transitions from a length on the order of hundreds to thousands of repeating unit cells (RUCs) to that on 
the order of a representative volume element (RVE) composed of enough microstructural elements (grains, 
inclusions, etc.), such that the RVE behavior is typical of the bulk composite (with an embedded localiza- 
tion), on average. Typically, an RUC does not have a physical length associated with it, but an RVE must 
have a physical length associated with it, and the RVE must contain a large enough volume that it captures 
the essence of the microstructure from a physical standpoint. Ref. 8,9 have studied the effect of the RVE 
size on deformation. However, when there is localization, the adequate RVE size will be different. Thus, 
the progression of strain localization is affected directly by the microscale fields in the RVE. For example, 
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if the RVE is in a compressive load environment, fiber kink-banding (localization) within the RVE leads to 
a post-peak softening response influenced by the size of the RVE as shown in studies reported by. 10 If the 
RVE is too small, a constraining effect leads to a kink band that shows a lower peak strength when compared 
to the asymptotic value of peak strength associated with a larger RVE that leads to the formation of an 
unconstrained and “free” kink band. 

Multiscale modeling is a popular technique for incorporating microscale effects into a structural scale 
model. With these methods, micromechanics models are linked to structural models and localization de- 
tails, at the microscale, are captured directly during down-scaling, while homogenization is employed for 
up-scaling, to communicate the effects of localization by smearing over the subscale model. The linking of 
scales can be achieved in a hierarchical, concurrent, or synergistic sense. 11 With hierarchical multiscale ap- 
proaches, micromechanics or subscale simulations are preformed a priori, and the results obtained from those 
simulations are utilized in subsequent macroscale, or structural level, models. With concurrent multiscale 
modeling, both the micro and macro scales operate simultaneously in time and space. Finally, synergis- 
tic multiscale models operate concurrently in time and hierarchically in both spatial scales (up-scaling and 
down-scaling), or vice versa. 

The micromechanics models employed can be either analytical, semi-analytical, or fully numerical. Typ- 
ically, analytical methods utilize mean field theories which are very efficient but only offer a single set 
of fields for each constituent. 12, 13 Fully numerical methods use the finite element method (FEM) or the 
boundary element method to model the microscale. 14-16 These methods can be computationally expen- 
sive; however, numerous authors have developed techniques to simplify the subscale FEM problem, yielding 
significant speed-up, for instance, we cite,. 7, 16-22 Semi-analytical methods offer the best balance between 
computational efficiency and solution fidelity. 23-29 There exist a plethora of multiscale techniques for fiber- 
reinforced composites in the literature, some of which, relevant to the present paper, have been reported in 
Ref. 11,29-35. 

Physically, a continuum material must posses a positive-definite tangent stiffness tensor, and, in fact, 
at a small enough scale, the material tangent stiffness tensor always remains positive-definite in order to 
satisfy the necessity of having a real, local speed of sound. 36 However for practical purposes, engineers must 
model structures at scales much larger than the characteristic flaws in the material. The homogenized (over 
a representative volume whose characteristic length is larger than the typical flaw size) continuum represen- 
tation of a material containing the nucleation and propagation of discontinuities, such as cracks or voids, 
will exhibit post-peak strain softening in the macroscopic stress-strain response. Although micromechanics 
and multiscale methods can be utilized to introduce details of the composite microstructure into structural 
analyses, erroneous numerical failure predictions can be obtained when post-peak strain softening is exhib- 
ited in the material, regardless of the scale. Loss of positive-definiteness of the tangent stiffness tensor leads 
to a material instability, which manifests as a localization of damage into the smallest length scale in the 
continuum problem. 36 Since the post-peak stress-strain relationship prescribes the energy density dissipated 
during the failure process, the total amount of energy dissipated is proportional to the size of the localization 
element, and in the limit as the element size is decreased, zero energy is required to fail the structure. 37,38 
Thus, the computational results will become pathologically dependent on the mesh size. 

A simple way to overcome this deficiency is to judiciously scale the post-peak softening slope by the 
characteristic element length such that the total energy release rate in the post-peak regime, after reaching 
a state of zero stress, is equal to the critical energy release rate, or fracture toughness, of the material (a 
fixed, experimentally obtained material parameter). This approach, known as the crack band or smeared 
crack approach, has been used by many authors to model strain localization due to failure within an FEM 
framework in a mesh objective manner. 36,39-44 Recently, the authors, Ref. 45,46, have implemented the crack 
band model within the generalized method of cells (GMC) 23 and high-fidelity generalized method of cells 
(HFGMC) 25 micromechanics theories and the predictions were verified against analogous FEM models. Non- 
local, or gradient-based theories have also been shown to prevent strain localization, eliminating dependence 
of the numerical solution on the size of the elements, as discussed in,. 47-49 However, the latter techniques 
require higher-order numerical interpolations and can be challenging to implement. 

Once a mesh objective theory has been implemented in a single scale analysis, the energy dissipation 
is preserved regardless of the size of the mesh, at that scale. However in a multiscale analysis, careful 
attention must be paid to how energy is being preserved and transferred across the various length scales. 
Ref. 50 commented that many multiscale models fail to define an appropriate localization limiter at the 
macroscale or fail to deliver any localization limiter to the macroscale from the microscale. If this is not 
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achieved energy dissipation is not necessarily preserved across the scales. However, if the appropriate theories 
are in employed at the correct scales and a consistent “handshaking” methodology is used that yields an 
appropriate localization limiter at the macroscale, then both the energy density and energy release rate will 
be preserved across the scales. 

The objective of the present paper, therefore, is to describe a method, using a two-scale analysis, that is 
energy preserving (in the sense of dissipation in the post-peak regime) at both scales with the requirement 
of consistent characteristic length relationships across the two scales. The method described is general, and 
its extension to multiple-scale (more than two length scales) analysis will become apparent. In Section IV, 
the effects of judiciously scaling the microscale RUC size to correspond to the macroscale element size are 
presented for two different multiscale example problems. 


II. Discretization Objective Progressive Failure Modeling at the Microscale 

The objective of multiscale modeling is to incorporate the influence of the material microstructure, at a 
suitable length scale, on the overall macroscopic response of the composite structure. To achieve this feasibly, 
an efficient microscale model is required. However, some fidelity at the microscale must be maintained if 
multiscale modeling is to realize any benefit. GMC offers an excellent balance between efficiency and fidelity, 
and is ideal for multiscale modeling. Since the method is semi-analytical, it offers computational speed that 
is superior to FEM based models (orders of magnitude faster). Additionally, GMC can be used to calculate 
local field distributions (essential for nonlinear problems involving inelasticity and damage), as opposed to 
mean field theories, which only offer a single field value (some chosen measure of the average) per constituent. 

Along with an appropriate microscale theory, physics-based constitutive laws (for the constituents in the 
composite) must be implemented for predictive capability. GMC admits a variety of linear and nonlinear 
constitutive laws, as well as any number of constituents. If the constitutive response of each constituent ex- 
hibits positive-definite behavior, then any suitable theory can be readily implemented within GMC. However, 
if the constitutive law exhibits post-peak strain softening behavior (i.e. , a negative tangent stiffness), then an 
appropriate formulation must be used to ensure objectivity with respect to refinement of the discretization 
used to represent the microstructure. 

One way to alleviate dependence on the level of discretization is to scale the energy density dissipated 
during the failure process, such that the energy release rate is preserved, independent of the size of the 
discretization element. One such model, the crack band theory, 39 was implemented within the GMC and 
HFGMC micromechanics theories to alleviate dependence of the failure solution on the size of the subcells 
used to discretize the RUC, as described earlier in Ref. 45,46. This implementation was utilized in this work 
to eliminate pathological dependence of the dissipated energy on the size of the subcells at the microscale. 


II. A. Implementation of the Crack Band Theory within GMC 

II. A. 1. Physical behavior of crack band 


The crack band model is intended to capture the behavior of a region of a material wherein numerous 
microcracks have initiated, and they eventually coalesce to form a larger crack. Figure 1 displays a crack 
band of width w c embedded in a continuum. The domain of the crack band is denoted as f Y and the 
remaining continuum as C l. The crack band is oriented within the continuum such that, for a given point 
within the crack band, the unit vector normal to the crack band is n. 

The total energy dissipated during the failure process is smeared over Cl', and the size w c of Cl' is a 
material property (characteristic length) directly related to the material fracture toughness. 39 


2 Qc 


1 

Ej 1 


-l 


(1) 


where ac is the critical stress for initiation of the post-peak regime in the ID material stress-strain law, 
and Et is the negative tangent slope in that regime. The fracture toughness Qc > or critical strain energy 
release rate, of the material is given by the area under the ID traction-separation law that governs the 
cohesive response of the separation of crack faces as a crack propagates in the material. The energy density 
dissipated during failure H-'V is related to the material fracture toughness through the characteristic length 
in the material. 

Qc = w c Wf (2) 
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II. A. 2. Formulation of crack band model 


Figure 2 shows a discretization of the continuum displayed in Figure 1. A magnified view of the crack band 
embedded in a single GMC subvolume (i.e., subcell) is also displayed in Figure 2. Since all of the energy 
dissipated in the crack band is smeared over the subcell volume, the subcell must be large enough to contain 
the crack band of width wc ■ Note that Figure 2 shows a 2D geometry for illustrative purposes, but the 
crack bands can also evolve in a general 3D setting. 

In GMC, a doubly-periodic RUC is represented with Np x iV 7 subcells. A GMC representation of a 
unidirectional, fiber-reinforced composite containing square-packed fibers is shown in Figure 3. Although a 
2D representation is shown, the calculated, local fields in the RUC are fully 3D, but there is no variation in 
the fields along the aq-direction. Subcells in the RUC are identified by the indices /3 (in the ^-direction) 
and 7 (in the aq-direction), and the RUC dimensions are given by H and L. 

The orientation of the crack band in subcell is given by the vector n^ 7 ' ) (see Figure 2) and is deter- 
mined from the local principal stress state In a monolithic material, it is postulated that 

cracks orient such that the crack tips are always subjected to pure mode I (opening mode) conditions unless 
there are constraints that restrict the crack orientation. These constraints can be the interface between two 
constituents, the interface between two adjacent plies, or a compressive loading state. Utilizing micromechan- 
ics, the constituents within the composite material are modeled discretely as separate, monolithic materials. 
Thus, in the pure matrix, there is nothing to constrain the crack band, and it can be assumed the crack band 
orients perpendicular to the principal stress with the largest magnitude, Id^ 7 ^ > > | d^ 7 ^ | , 

if > 0 (tensile). Under these conditions, a crack oriented as such, is subjected to locally pure mode I 
loading, which is the most energetically favorable state. Although, the resulting global behavior may appear 
to be mixed mode because of the influence of the fibers on the local stress state driving the matrix crack 
band path. Crack band initiation is determined using a very simple, but physically-based criterion . 3,4 


df 7) 

(dr) 

a C 


= 1 , 


A/U) 


' l ^ 


> 0 


( 3 ) 


where is the cohesive strength of the matrix material. Once the crack band has initiated, the crack band 
orientation is assumed to be fixed within a particular subcell, as time evolves. If the direction of maximum 
principal stress is deemed an unsuitable choice for the crack band normal, other directional measures can be 
readily employed to determine the crack band orientation. 

Once the orientation of the crack band has been determined, the subcell compliance is rotated into the 
principal frame using the transformation matrix. 

T =[n^ 7) n^ 7) n^ 7) ][e ie 2 e 3 ] (4) 

where n^ 7 \ n^ 7 \ and are the principal stress directions, and e±, e 2 , and e 3 are the unit basis 

vectors. All material degradation due to crack band evolution is imposed on the rotated compliance 
the components of which are given by: 


q(A y) — t .t . T, fAl 

°ijkl ± pi- L q3° pqrs ±kr- L ls ) 

The strain energy released during the formation of new surfaces corresponding to the growth of cracks 
within the crack band is assumed to be dissipated over the entire subcell volume. Therefore, the post-peak 
softening slope Ej ^ , and the strain at which the principal stress state is zero , are calculated using the 
characteristic length of the subcell /[A and the material fracture toughness 
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where U 110 is the undamaged, axial Young’s modulus in the principal frame. The characteristic length of 
the subcell is determined as the dimension of the subcell running parallel to (see Figure 2). 
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It should be noted that e\ ^ r ' > must be less than zero; therefore, by Equations (6) and (7), a restriction 
is placed on the maximum allowable subcell size. 


(07) 


< 


2 

,A07) 2 


( 8 ) 


Use of a subcell violating Equation (8) results in non-physical snap-back. If this is unavoidable, the strength 
of the material can be scaled to accommodate size effects, and a sudden drop in the stress to zero can be 
employed subsequent to failure initiation, as descried in,. 39 In addition to a maximum size restriction, the 
subcells must also be larger than the actual characteristic length of the material wc- 


WC < # 7) 


(9) 


If the subcell size is smaller than wc, then non-local theories must be used because it is unrealistic for the 
failure to localize to a scale smaller than the characteristic length of the material; thus, localization to a 
single subcell must be prevented outright. 


The local, rotated, subcell strain state e\ 
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is used to degrade the rotated compliance components. The scalar, mode I, damage factor D^ 3 ' 1 ' 1 is calculated 


using the rotated strain corresponding to 
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where e^ 5 * 7 ' 1 is the value of when the initiation criterion, Equation (3), is satisfied. If Dj 3l ' ) is less than 
zero, no damage occurs, while a maximum damage level of one corresponds to a zero stress state on the 
softening stress-strain curve. Also, damage healing is inadmissable. 
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> 0 
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The components of the rotated compliance matrix are degraded with the damage factor according to: 
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Initially, the crack band is free of any shear tractions. To prevent numerical instabilities, the crack band 
orientation is fixed for a particular subcell, upon initiation. Thus, the s[ g 7 g and <§(212 shear compliances in 
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the rotated frame are degraded, as well as the fqf)/ compliance, so that the faces of the cracks within the 
band normal to n ^ 3 ~ 1 ^ are free of normal and shear tractions when all of the crack band energy has been 
dissipated (i.e. = Qj^)- scales with Z^ 6 7 ' 1 in a discretized continuum. In reality Wf in 

Equation (2) is fixed; however Wp jy \ which is smeared over the subcell volume, changes as a function of 
the characteristic subcell length. A mixed-mode law could be introduced to ensure that as shear tractions 
develop, the appropriate mode II strain energy release rate is dissipated. However, most mixed-mode theories 
utilize an initial mode mixity parameter to calculate the effective traction-separation law. 51,52 For a crack 
band initially oriented perpendicular to the maximum principal stress direction, this mode mixity parameter 
would be zero. 

Once the compliance in the rotated frame is degraded, it is transformed back to the global frame to yield 
the new compliance of subcell (/Uy). 


o(fil) rp—lrji—1 a(0j)rp—lrp—l / 1 , \ 

J ijkl ± pi ± qj J pqrs ± hr ± Is GU 

Although only one direction is damaged in the principal frame, when transformed back to the global frame 
any general direction may accumulate damage. Additionally, damage introduced in the principal frame, 
through Equation (13), can induce normal-shear coupling in the global frame. 

If the maximum principal stress is compressive, then any initiated crack band closes and the rotated 
compliance tensor is modified to reflect this. 
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(15) 


Damage can no longer evolve under compression conditions, and the stiffness in the direction perpendicular 
to the crack band reverts to its undamaged value. The shear stiffness in the plane of the crack band remains 
at the previously attained, maximum degradation state. Ref. 45,46 implemented a Mohr-Coulomb initiation 
criterion, and mode II evolution law, to allow shear cracks to develop under a principally compressive stress 
state. This formulation is omitted from this paper for brevity. The focus of this paper is to demonstrate the 
necessity for energy preservation across the scales. An example is chosen that ensures local, tensile stress 
states in the principal directions, and no prediction of compressive failure is attempted. 


III. Multiscale Modeling of Composite Structures Using GMC 

RUCs modeled with GMC can be easily linked to higher-level, structural FEM models. The integration 
point strains, from the FEM model, are applied to the RUC and the local subcell fields are determined 
using GMC; this process is referred to herein as clown-scaling for clarity, although in literature it is more 
often referred to as localization. If the subcell material behavior is nonlinear, the local stresses and strains 
are used to calculate the local stiffnesses, inelastic strains, thermal strains, and/or state variables. Since 
the local stresses and strains depend on the local stiffnesses, inelastic strains, and thermal strains, some 
iterations of this procedure may be necessary to resolve the correct local fields if there is a high degree of 
nonlinearity. The RUC is then homogenized and the global stiffnesses, inelastic strains, thermal strains, 
and/or state variables are computed and passed onto the macroscale in a process referred to as up-scaling. 
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The global stresses and material Jacobian at the integration point are then computed using the up-scaled, 
global, homogenized fields. A schematic showing this procedure and the corresponding scales is presented in 
Figure 4. 

Since failure is dictated by extreme values, it is hypothesized that the fiber-matrix architecture influ- 
ences the evolution of failure, as it influences the local extreme values. The extreme values are a result of 
non-homogeneity due to property mismatch between constituents, and geometrical packing details of the 
constituents. Therefore, failure is modeled at the microscale, within GMC, using the crack band approach 
presented in Section II. Homogenization of the microscale is used to influence failure mechanisms at the 
higher scales (in this two-scale case, the higher scale being macroscale). 


III. A. Preserving Energy Across Scales Using Consistent ‘Handshaking’ Methods 

In order to retain mesh objectivity in a multiscale model, both the energy density and energy release rate 
must be preserved at, and across, all scales. Utilizing the crack band model (see Section II) within the GMC 
micromechanics ensures objectivity at the microscale. However, there must be a consistent relationship 
between the characteristic lengths at each scale in order to preserve the energy release rate. 

Communication between the finite element integration points and the GMC RUC is achieved through 
the global stress a and strain e measures, see Figure 5. The global stress and strains are the fields at 
corresponding element integration points and also represent the average global stresses and strains applied 
to the RUC. In addition the material Jacobian, or tangent stiffness matrix, of the RUC calculated by GMC 
can be passed to the FEM and used to formulate the element stiffness matrix and/or estimate the next 
global strain increment. 

The energy density can be defined at the microscale for the RUC 


W RUC 


^RUC j-RUC 
a ij de ij 


(16) 


where <J RUC and e RUC are the components of the global, average stress and strain tensors applied to the 
R.UC, and the energy density can be defined at the macroscale for an integration point 


W int 


-int 

’ ij Utij 


(17) 


where ct“ 4 and e)” 4 are the components of the stress and strain tensors at the finite element integration 
point. 

Since, using the previously described handshaking method, the global stresses and strains are maintained 
across the scales (i.e., dj” 4 = & RUC and e)" 4 = e RUC , it follows that the energy density is automatically 
preserved. 

W int = W RUC (18) 

The energy density for the finite element can be taken as the volume averaged sum of the energy density of 
all the integration points within that finite element 

1 nint 

w e = — - ^2 W int V int (19) 

int — 1 


where V mt and V e are the integration point subvolumes and element volume, respectively. 

In addition to the energy density, the energy release rate must be preserved across the scales. The energy 
release rate, or fracture energy, is defined as the energy density (in the absence of plastic strain accumulation) 
multiplied by some characteristic length associated with the localization band. 36, 39,53-55 Thus, the total 
energy release rate for the RUC at the microscale is given by 

QRUC = tRUCyyRUC ( 20 ) 

where l RUC is the characteristic material length associated with strain softening in the RUC. Similarly, the 
energy release rate for the integration point at the macroscale is defined as 


G 


int 


= i i g t w int 


(21) 
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where Z™* is the characteristic length associated with the integration point. It follows that the total energy 
release rate of the finite element is obtained from the volume averaged sum of the energy release rates of the 
integration points. 

nint 

g e = — ^2 g int v int ( 22 ) 

int=l 

In order for the energy release rates at both scales to be equal, the characteristic lengths associated 
with the finite element integration point and the RUC must equivalent. At the subcell level, it is trivial to 
define the characteristic length, (see Section II. A) — it is simply the length of the subcell perpendicular 
to the crack band. However, the characteristic length of the RUC is a result of the evolution of the crack 
band(s) in multiple subcells contained within the RUC. Since, the crack band path is complex, and the final 
configuration is not known a priori it cannot be calculated and delivered to the macroscale to ensure that 
the lengths are equal. 

Moreover, there is a restriction placed on the minimum allowable size of the macroscale element. The 
element (or integration point subvolumes) must be large enough to contain the characteristic length of the 
R.UC, or a suitably large RVE containing enough fibers needed to represent the localization objectively. 8, 9,46 
If the element is too small, then microscale homogenization techniques are not valid and physics is violated. 

It is hypothesized here that if the size and shape of the RUC and integration point subvolume are 
identical, then all lengths must be consistent across the scales, automatically. This statement holds true for 
the characteristic lengths also, even if they cannot be calculated or defined. Unfortunately, GMC only admits 
rectangular RUCs (although a new isoparametric mapping technique has been developed for HFGMC 56 ), and 
it is impractical to assume the finite element mesh will contain only rectangular elements. Particularly when 
modeling structural components, generally shaped quadrilaterals for 2D, or hexahedra for 3D, structures are 
more appropriate . 

It can be assumed that the characteristic length of a 2D element integration point is equal to the square 
root of its area. 

lint _ yjAint (23) 

A similar relationship can be assumed for the doubly-periodic RUC 

iRUC = JJruc (24) 


where H and L are the RUC dimensions given in Figure 3. If the dimensions H and L are scaled such that 
the new values (denoted by a prime symbol) are given by 


H' = 


L' = 


l^H 

Vhl 

(25) 

l^L 

Vhl 

(26) 


then utilizing Equations (25) and (26) in Equation (24) results in equivalent characteristic lengths at both 
scales, that is, 

1% UC = Zj?* (27) 

Finally, substituting Equations (18) and (27) into Equations (20) and (21) yields 


G 


int 


gRUC 


(28) 


which, along with Equation (18), is a necessary condition for mesh objective multiscale modeling of strain 
softening behavior. 


III.B. FEAMAC Multiscale Framework 

In Section IV, a multiscale framework is utilized to perform simulations of fiber-reinforced composite struc- 
tures by modeling the fiber-matrix architecture as an RUC at the microscale using GMC and linking the 
microscale to the lamina/laminate level (macrsocale) FEM model. A synergistic approach is employed that 
executes concurrent multiscaling in time, but two-way hierarchical multiscale in space. 11 The commercial 
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finite element software, Abaqus 6.11-1 57 is used as the FEM platform, and the MAC/GMC core microme- 
chanics software 58, 59 is used to perform microscale calculations, and the scales are linked using the FEAMAC 
software implementation. 27 The number of subcells and constituents within an RUC at the microscale is 
completely general and as many, or as few, may be used as needed to accurately characterize the micro- 
architecture of the composite. Moreover, any constitutive model can be used in the subcells of the RUC to 
represent the mechanical or thermal response of the constituents. 

FEAMAC consists of four Abaqus/Standard user defined subroutines, 57 as well as six subroutines exclu- 
sive to the FEAMAC package (see Figure 6). The Abaqus/Standard user material UMAT subroutine provides 
the strains, strain increments, and current values of state variables to MAC/GMC through the front end 
subroutine FEAMAC. For this work, crack band calculations are performed within MAC/GMC, and FEAMAC is 
called iteratively until convergence is acheived at the microscale. MAC/GMC then returns a new stiffness and 
stress state, as well as updated state variables, to the UMAT via the FEAMAC subroutine. The Abaqus/Standard 
user subroutine SDVINI initializes the state variables used in UMAT. The Abaqus/Standard user subroutine 
UEXPAN is used for thermal analysis by providing the integration point temperature, temperature increment, 
and current state to MAC/GMC, which in turn, calculates new thermal strains and thermal strain rates. 
Problem set-up task, initialization, and writing MAC/GMC level output data to files is achieved through 
the Abaqus/Standard user subroutine UEXTERNALDB, which evokes communication among Abaqus/Standard 
and the FEAMAChPRE and FEAMAC_PL0TS subroutines. The reader is referred to Ref. 27 for further details 
on the FEAMAC software implementation. Furthermore, a wrapper (FEAMAC/Explicit) was developed to 
link FEAMAC to the Abaqus/Explicit FEM software through a VUMAT Abaqus user material subroutine. 60 

IV. Numerical Examples 

IV. A. Example 1 - Finite-notched DCB 

IV.A.l. Multsicale Model Details 

To demonstrate the effects of appropriate microscale RUC length scaling, a simple multiscale example was 
chosen representing a modified DCB specimen. The global geometry of the modified DCB is shown in Figure 
7. The total specimen length was 100 mm, with a height of 6 mm. Plane strain conditions were assumed. 
A 25.5 nun notch was placed at one end of the specimen. Since continuum elements were used to model 
the domain and incorporate damage, a notch with a tip of a finite diameter (1 mm) was chosen so that a 
converged stress solution could be obtained. If a sharp crack had been chosen, the stress concentration at 
the notch tip could not be calculated properly. The focus of this work was to isolate mesh dependence due 
to post-peak strain softening. Hence, it was desirable to eliminate the dependence of the elastic solution on 
the global mesh size, altogether. 

The entire domain was meshed using 2D, plane strain, reduced integration, quadrilateral CPE4R, and 
triangular CPE3 Abaqus elements. 57 A mesh sensitivity study was conducted and it was determined that 
using an element size l e (see Figure 7) of 0.075 mm, or smaller, ahead of the notch tip would yield a converged, 
elastic, stress state at the notch tip. Figure 8 shows the coarsest mesh used in the studies incorporating a 
0.075 mm by 0.075 nun elements ahead of the notch tip. Figure 8b shows a magnified view of the mesh 
surrounding the notch tip. Assuming a fiber diameter of 5 /jm, an element of this size would contain nearly 
170 fibers, which is adequate to justify the use of periodic boundary conditions, initially. 

The macroscale domain lies in an x-y-z coordinate system. The model is intended to simulate matrix 
cracking in a 90°, IM7-8552 laminate, so the fiber direction is aligned with the global z-axis. For simplicity, 
and efficiency, the majority of the elements, represented with green in Figures 7 and 8, utilized a single-scale, 
elastic, transversely isotropic constitutive law. The elastic properties for IM7/8552 are given in Table l. 43 

The elements along the expected macroscale crack path, represented in yellow, were linked to a microscale, 
GMC RUC using FEAMAC. The microscale model consisted of a 7 subcell by 7 subcell RUC with a local X\- 
X 2 -X 3 coordinate system. The RUC was composed of 13 subcells comprised of the fiber constituent (colored 
blue, see Figure 7) and 36 subcells comprised of the matrix constituent (colored green). A fiber volume 
fraction Vf of 60% was maintained. The axial, fiber direction, x\, is aligned with the global z-axis. This 
choice of RUC may not be completely suitable once localization occurs, since the size of the localization 
typically spans one, or multiple, unit cells. An RVE containing multiple fibers is more appropriate to 
capture the effects of the microscale localization objectively. The intention of this work was not to offer 
blind failure predictions using the multiscale methodology, or to investigate solutions objectivity relative to 
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RVE size and complexity, but rather to demonstrate the necessity for energy preservation and consistent 
handshaking methods across the scales. Thus, the simple and efficient RUC was chosen for its computational 
speed advantages and tractability, so that numerous multiscale simulations could be achieved in a reasonable 
amount of time. Furthermore, determining an adequate size and level of detail for the RVE, when there is 
localization at the subscale, remains an ongoing research area. 8, 9,46 

The elastic properties of the fiber (transversely isotropic) and matrix (isotropic) are given in Table 2. 
These properties were calibrated such that the global RUC properties calculated using the 7x7 GMC RUC 
corresponded to the elastic properties given in Table 1. The transverse fiber stiffness E 22 , and axial shear 
modulus G\ 2 were taken from Ref. 61. 

The crack band theory, detailed in Section II. A, was used to dictate failure in the matrix subcells. 
The transverse strength and mode I fracture toughness of a unidirectional, IM7/8552 ply was reported 
by Ref. 43 as 62.3 MPa and 0.2774 kJ/m 2 . The same, single-fiber, 7x7 RUC was used, assuming a fiber 
diameter of 5 /iin, to calibrate the matrix strength and fracture toughness. Using the standalone MAC/GMC 
micromechanics software, The R.UC was loaded under transverse strain and the matrix properties were 
adjusted until the global RUC exhibited a strength and total energy release rate upon failure corresponded 
to the values reported in Ref. 43. The resulting properties are presented in Table 3, and the stress strain 
response of the R.UC is plotted in Figure 9. Note that the incredibly high strain to failure (> 1.5) is a 
result of the crack band opening displacement smeared over the subcell, not elastic strain in the material 
(continuum). Since the single fiber R.UC utilized realistic fiber dimensions, a reasonable opening displacement 
yields a seemingly very large homogenized strain. At larger, structural scales, the same crack band opening 
displacement would produce more reasonable homogenized strains. 

The focus of this work was to determine the effects of incorporating a consistent length across the scales in 
a multiscale analysis. The yellow domain containing multiscale elements, shown in Figures 7 and 8, utilized 
fixed element shapes and sizes. Since single integration point elements were used, the characteristic length 
of the elements and integration point subvolumes are identical. 

l e c = ( 29 ) 

The element sizes in front of the notch tip corresponding to the multiscale elements were decreased from 
0.075 mm, the maximum allowable element size needed to obtain a converged stress state at the notch tip. 
Two types of analyses were performed: scaled and unsealed. Scaled analyses utilized Equations (25) and (26) 
to scale the dimensions of the RUC so that they corresponded to the element dimensions. Unsealed analysis 
assumed fixed R.UC dimensions equal to 0.075 mm x 0.075 mm. When the RUC dimensions are scaled the 
volume of the RUC represents a physical volume, similar to an RVE. However, due to the simplicity of the 
fiber/matrix architecture it cannot statistically represent the heterogeneous nature of the composite and 
periodic boundary conditions are still utilized; thus, this scaled unit cell will still be referred to as an R.UC. 

In addition to two types of analyses, two types of elements were used along the notch tip front. Square 
elements were used as a baseline because the assumed characteristic length Iq (Equation (23)) and the 
dimensions of the element sides l e (see Figure 7) are equivalent. Therefore, if RUC length scaling is performed, 
the dimensions of the R.UC and the element are equal. Thus, every length definition is maintained across 
the scales and Equation (27) automatically holds without needing to make the assumption in Equation (23). 
With a consistent size and shape at both scales, the characteristic lengths do not need to be defined explicitly, 
yet it can be accepted that the appropriate localization limiter will be delivered to the macroscale (due to 
the equivalence of all length definitions) and Equation (28) will be satisfied assuredly. 

In addition to square elements, triangular elements are used (see Figure 8c) wherein Iq ^ l e . Therefore, 
the exact size and shape were not preserved across the scales, only the volume. Hence, not all length 
definitions are inherently preserved. The results from these scaled simulations will indicate if the assumption 
made with Equation (23) is always valid. 

IV. A. 2. Numerical Results 

Figure 10 shows the ultimate load predicted from the two different analysis types (scaled and unsealed) 
utilizing two different element shapes (square and triangular). The element size was varied from 0.075 mm 
to 0.045 mm (which could still accommodate 61 fibers at a 60% fiber volume fraction) in increments of 
0.05 mm. Results from analyses using square element are indicated with square symbols, and results from 
analyses using triangular elements are indicated with triangular symbols. Filled markers signify that the 
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RUC dimensions were scaled according to the assumed characteristic element length (Equation (23)) using 
Equations (25) and (26), and open markers denote that a fixed RUC size of 0.075mm x 0.075mm was used. 

When square elements with scaled RUCs were used, the ultimate load predicted during the simulations 
remained within 1% of the mean 8.02 N. Whereas, not performing the scaling appropriately led to a steady 
reduction in the ultimate load as the element size was decreased. Using a square element with length 0.045 
mm, and an unsealed RUC, yielded an ultimate load prediction of 7.02 N, a 13.5% error from the mean 
of the scaled analyses. This pathological mesh dependence is the most significant limitation of multiscale 
modeling identified by Ref. 50, as it yields inconsistent and unreliable predictions. However, as shown, 
with the appropriate “handshaking” methodology and the use mesh independent damage theories at the 
microscale, this pathological dependence on the discretization size can be altogether eliminated at all scales. 

Scaling also eliminates pathological mesh dependence when using triangular elements as shown in Figure 
10. The ultimate load predictions vary from the mean by 1%, at most. However, the mean ultimate load 
obtained from the simulations utilizing triangular elements varies from the mean ultimate load results calcu- 
lated using square elements by 9.4%. This significant error indicates that the estimate of the characteristic 
length using Equation (23) is incorrect for triangular elements. When the scaling is not implemented, the 
ultimate load prediction pathologically decreased as the size of the triangular elements was reduced resulting 
in a maximum error from the mean of 21% with an element size of 0.045 mm. 

It should be noted that in Figure 10, unsealed results obtained using both element types are comparable. 
Since a fixed RUC size is used, the effects of the scaling are eliminated from the problem, and the displayed 
trends are purely a result of the change in the element size. This indicates that a better estimate of the 
localization limiter may have been the width of the advancing macroscale softening band (localization), 

l e c = l e (30) 

not the assumed characteristic element length computed using Equations (23) and (29). 

The ultimate load predictions obtained when the RUC dimensions are scaled according to Equation (30) 
are displayed in Figure 11 along with the previous results obtained from the scaled simulations (utilizing 
both square and triangular elements). The lines represent the mean value for the corresponding simulations. 
It can be observed that using the element length rather than the square root of the area of the element 
improves the results for the simulations utilizing triangular elements. The mean error using square elements 
improved from 9.4% to 1%; however the results are not identical. 

For the problem demonstrated, where the macroscale localization growth is self-similar, it may seem that 
a suitable localization limiter would be the width of the macroscale localization band. However, if the crack 
band path at the microscale is not known beforehand it is impossible to assume the localization limiter at 
the macroscale because it is inherently coupled to the localization length(s) at the microscale. 

IV. B. Example 2 - 2D Plane Strain Square Section 

IV.B.l. Multsicale Model Details 

A second example, first presented in Ref. 62, is included in this manuscript to further demonstrate the 
necessity for appropriate “handshaking” in order to yield a mesh objective and energetically accurate solution. 
In this multiscale example, the macroscale was composed of a square domain meshed using plane strain 
quadrilateral elements (see Figure 12) or plane strain triangular elements (see Figure 13). As in the previous 
example, square elements are used to establish a baseline solution wherein size and shape of the finite 
elements and microscale RUCs are identical; thus, preserving all other possible definitions of length across 
the scales. 

Failure is introduced into the matrix subcells of the microscale RUC, and the effects are transferred to the 
macroscale through an overall reduction of the stiffness of the RUC (due to a local reduction in the stiffness 
of the failing matrix subcells). The crack band model presented in Section II is used to govern failure at the 
microscale. Only the yellow and red elements in Figures 12 and 13 were modeled with multiscale FEAMAC 
elements. The green elements were modeled using the Abaqus elastic material model. Free edge effects and 
the accuracy of the local fields surrounding the failure path front are influenced by the mesh refinement, 
to some degree, which may affect the failure path in the simulation. Therefore, failure was restricted to 
the yellow and red elements to simplify the problem as much as possible and try to ensure that the energy 
dissipated in the model is a result of the multiscale methodology, not other influences. Moreover, the matrix 
subcells of the RUC located in the red element were given a slightly lower cohesive strength than the matrix 
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of the RUCs within the yellow elements. This was imposed to imitate an imperfection and facilitate the 
localization of initial failure into a single element. Without this imperfection, all of the FEAMAC elements 
would fail simultaneously because the stress and strain fields are uniform throughout the FEM mesh. After 
failure initiates in the RUC located in the red element, failure progresses into the RUCs located within the 
adjacent yellow elements until the failure reaches the free edge of the FEM mesh. 

The microscale RUCs are linked to the macroscale integration points such that the global (macro) z- and 
local (micro) aq-axes are aligned. Thus, the fibers can be considered running “out of the page.” The local, 
transverse, ^-direction is coincident with the global x-direction, as is the the local X 3 with the global y. 

Reduced integration elements were used in both the square and triangular element cases. Therefore, 
Equation (29) holds. Furthermore, the definition of characteristic length defined by Ref. 57 is utilized 
(Equation (23)). Two types of multiscale simulations were performed: scaled and unsealed. For the scaled 
analysis the dimensions of the RUC were scaled according to Equation (27), assuming Equation (24). Un- 
sealed simulations incorporated a fixed length for the sides of the RUC equal to 5.98 /tm. 

Several simulations were performed with various levels of macroscale mesh refinement. In all simulations, 
the left edge of the macrsocale domain was fixed in the ^-direction, while the bottom edge was fixed in the 
//-direction, and a uniform displacement was applied to the right edge in the x-direction. The domains were 
loaded until final failure; i.e. , a complete loss in load carrying capability. 

IV.B.2. Numerical Results 

Figure 14 shows the strain energy release rate ( Q ) as a function of macroscale characteristic element length at 
final failure for the numerous different simulations. Simulations performed using square elements are shown 
with square data points, while triangular elements are indicated with triangular symbols. Filled symbols 
denote that a scaled analysis was performed, while unfilled symbols mark unsealed analyses. 

The element length that is closest to the realistic size of the RUC (L = 5.98 /tm) is an element length 
of 5.99 n m (167 elements x 167 elements), and is only 0.17% larger than the actual RUC length. Thus, the 
scaled solution obtained with = 5.99 /tm is considered to be the most accurate solution. For this mesh, Q 
was calculated to be 0.227 kJ/m2; this value will be used as the baseline energy release rate for comparison 
with the other simulations. The most erroneous Q for the scaled simulations utilizing square elements was 
calculated from the solution using the coarsest mesh (32.3 /mi), and was 16% higher than the baseline 
solution; although, the element length was 438.7% larger. As the mesh was refined, Q approached 0.227 
kJ/m2. 

If the dimensions of the RUC were not scaled, the error observed in the energy release rate was very 
significant. The coarsest mesh using square elements exhibited 209.6% error from the baseline. Although, 
when the characteristic length of the integration point domain and the size of the RUC are comparable (i.e., 
5.99 /tm), Q only exhibited 1.6% difference from the baseline. However, the energy release rate calculated 
for the unsealed simulations dropped suddenly when the element size approached the size of the RUC. For 
example, the average Q for the next three largest meshes (6.14 /tm, 7.41 /tm, and 8.70 /tm) is 44.5% higher 
than the baseline. 

This example also demonstrates that scaling the size of the RUC such that it corresponds to the dimen- 
sions of the finite element minimizes the error. The variation that is observed in simulations with scaling, 
as the mesh is refined, can be attributed to local effects in the macroscopic FEM model such as an increase 
in the stress concentration ahead of the failure path front, as the mesh is refined. 

The simulations incorporating triangular meshes exhibited an overall larger degree of error (when com- 
pared to the baseline) in the strain energy release rate. This is to be expected, since the shape the element 
and the shape of the RUC are disparate. The characteristic element length, calculated using Equations (23) 
and (29), that was closest to the assumed actual RUC length used in the unsealed simulations, is 5.94 /tm. 
When the RUC dimensions were scaled, Q was calculated to be 21.6% higher than the the baseline calculated 
using square elements. The largest Q calculated (when RUC scaling was enforced), corresponding to the 
coarsest mesh with (33.7 /tm), was 74.9 % higher than the baseline. 

For the unsealed simulations utilizing a triangular mesh, the largest Q calculated corresponded to the 
coarsest mesh (33.7 /tm), and was 281% larger than the baseline. The smallest Q calculated was 23.8%, 
corresponding to an elements size of 6.15 /tm. The mismatch between the element shape and the RUC shape 
introduces some error when compared to the baseline square solution, just as with Example 1. Thus, the 
use of triangular elements in an FEAMAC simulation should be avoided. However, the error introduced 
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through refining the mesh can be minimized by scaling the RUC dimension by the characteristic length of 
the element. 


V. Conclusions 

A mesh objective crack band theory for modeling post-peak strain softening was implemented within 
the GMC micromechanics theory. Use of the proposed methodology is shown to lead to mesh objectivity in 
the macroscale model in a multiscale computation. In multisclae computations, there is no guarantee that 
mesh-objectivity at the subscale will automatically ensure mesh objectivity at the higher scales. To ensure 
mesh objectivity throughout, both the energy density and energy release rate must be preserved across all 
scales. To achieve this, some consistent, localization limiter must be delivered to the macroscale, as discussed 
in Ref. 50. 

The multiscale method proposed here uses the FEAMAC framework to link the MAC/GMC microme- 
chanics software (which contains implementations of GMC) to the Abaqus finite element software. In this 
synergistic methodology, integration point strains are passed to the microscale in the down-scaling step. 
Micromechanics calculations are performed, and the global RUC stiffness (material Jacobian) and stresses 
are passed back to the macroscale integration point in the up-scaling step. Equilibrium is ensured at each 
scale through iterations. 

It was hypothesized that if the dimensions of the microscale RUC were scaled according to a characteristic 
macroscale element length, then the ensuing computational results would be insensitive to changes in the 
macroscale mesh. It was assumed that the characteristic element length was equal to the square root of the 
element area (for 2D elements). This would ensure that, if square elements were used, both the size and 
shape of the R.UC and element would be identical. Thus, any possible length definition would be preserved 
across the scales. 

To test this hypothesis, two simple numerical examples were devised, and two types of analyses were 
performed: scaled and unsealed. Scaled analyses utilized the scaling procedure described previously, in lieu 
of fixed RUC dimensions; whereas unsealed analyses utilized the baseline RUC dimensions. Additionally, 
both square and triangular elements were utilized in the simulations, separately. If the scaling was enacted, 
the ultimate load predictions were insensitive to refinements in the macroscale mesh. However, if the RUC 
dimensions were fixed, the predicted ultimate load continued to drop as the element size was decreased, 
displaying pathological mesh dependence. This trend was true for both the square and triangular elements, 
but the scaled triangular elements produced a mean ultimate load higher than the square elements. This 
indicates that the assumed localization limiter used to scale the RUC for the triangular elements was de- 
pendent on the localization; therefore, the width of the global localization band was used to scale the RUC 
in the case of triangular elements. This gave improved results. However, the macroscale and microscale 
localization paths are heavily coupled and typically not known prior to simulation or testing of more com- 
plex structural components. Furthermore, the choice of an R.UC may not be appropriate to represent the 
microstructural behavior when localization occurs. Instead, a sufficiently large RVE may need to be used. 
Further studies must be conducted to determine the size and complexity of the RVE needed to capture 
microscale deformation localization objectively. 

The results of this study further support the argument that a suitable localization limiter must be 
delivered to the macroscale. 50 However, if the size and shape of the macroscale problem correspond directly 
to the size and shape of the macroscale element, then the correct localization limiter is maintained even it 
is not defined. Unfortunately, most realistic structural components cannot be meshed using rectangular or 
square elements. Therefore, the challenge remains to calculate a physically justified, fully general, macroscale 
localization limiter to be used to scale the R.UC. This is not trivial as the macroscale localization limiter 
depends inherently on the microscale localization path, which is not known in advance. 
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Property 

Value 

E zz 

171.4 GPa 

E 

J-'XX 

9.08 GPa 

Vzx 

0.32 

G zx 

5.29 GPa 


Table 1: Elastic transversely isotropic properties for IM7/8552 lamina used in single scale elements . 43 


8552 Matrix Properties 

Value 

IM7 Fiber Properties 

Value 

E m (calibrated) 

4.97 GPa 

E{ x (calibrated) 

286 GPa 

v m (calibrated) 

0.36 

fJ 61 
-^22 

12.4 GPa 



v{ 2 (calibrated) 

0.29 



^23 (calibrated) 

0.29 



ryf 61 
Lt 12 

20.0 GPa 


Table 2: Elastic properties of IM7 carbon fiber and 8552 epoxy matrix constituents used in microscale GMC 
RUC. 


Mode I Cohesive Strength 7) 

Mode I Fracture Toughness 

56.5 MPa 

0.266 kJ/m 2 


Table 3: Crack band failure parameters used in 8552 matrix subcells. 
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Q 



Figure 1: 


Crack band domain fl' of width w c oriented normal to vector n within a continuum Q. 



Figure 2: Crack band embedded in discretized continuum. Magnified discretization element displays crack 
band orientation, as well as, characteristic length of discretization element. 
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Figure 3: Doubly-periodic representation of a unidirectional, fiber-reinforced composite containing square- 
packed fibers using GMC. 



| HOMOGENIZATION 


GMC 


FEM 


LOCALIZATION 


’ dAs 


Figure 4: Schematic representation of the up-scaling and down-scaling steps employed in the synergistic 

multiscale technique utilizing GMC and FEM. 



Figure 5: Schematic communication between macroscale finite element integration point and microscale 

GMC RUG. Global strains are passed down to the microscale and global stresses are passed back up to the 
macroscale. 
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Figure 6: FEAMAC multiscale framework software implementation architecture. 



Figure 7: Multiscale model of finite-notched DCB specimen composed of a 90° laminate. The macroscale 
domain was modeled using traditional plane strain elements and the global x-y-z coordinate system. Elements 
within the green domain utilized a single-scale, transversely isotropic, constitutive law (^-axis represents the 
fiber direction). Elements within the yellow domain, ahead of the notch tip, were linked to a microcale GMC 
RUC. The microscale RUC consisted of a doubly-periodic, 7 subcell by 7 subcell RUC with a local X 1 -X 2 -X 3 
coordinate system. The RUC contained 13 fiber subcells (colored blue) and 36 matrix subcells (colored 
green), and the local fiber direction, x,\ . was aligned with the global 2 -axis. 
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(a) Global view. 



(b) Magnified view of square element mesh 
near notch tip. 



(c) Magnified view of triangular el- 
ement mesh near notch tip. 


Figure 8: FEM mesh of multiscale modified DCB. 



Figure 9: Global Transverse stress versus transverse strain response of RUC assuming a 5 /jm diameter fiber 
and 60% fiber volume fraction. Matrix failure parameters were calibrated so that global RUC strength and 
energy release rate upon failure equaled values reported in Ref. 43 for IM7/8552. 
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Figure 10: Ulitmate load predicted with multiscale, modified DCB simulations. Results compared for 

microscale RUC scaled by element length to fixed RUC dimensions of 0.075 mm x 0.075 mm. Results using 
square and triangular elements also compared. 
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Figure 11: Ulitmate load predicted with multiscale, modified DCB simulations. Results for microscale RUC 
scaled by assuming two different measures for characteristic element length. Results using square (previous 
results, from Figure 10) and triangular elements also compared. Lines represent mean values of corresponding 
simulations 



t 

Figure 12: Multiscale example problem. Macroscale FEM model consists of square, 2-D, plane strain, 

reduced integration CPE4R elements with dimensions le x le and global x-y-z coordinate system. The 
microscale model consists of a doubly-periodic, 7 subcells x 7 subcells GMC RUC with local xl-x2-x3 
coordinate system. The RUC contains 13 fiber subcells (colored blue), and 36 matrix subcells (colored 
green). The global dimension of the RUC are L x L, and the local subcell dimensions are h? x 1?. The 
red and yellow FEM elements are linked to the microscale RUC. The subcells within the RUC in the red 
element is given a slightly lower cohesive strength than those in yellow elements. The green elements in 
the macroscale FEM mesh are modeled as a standard, transversely isotropic, elastic material with elastic 
properties equal to the homogenized elastic properties of the GMC RUC. 
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Figure 13: Multiscale simulation utilizing triangular elements at the macroscale. 
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Figure 14: Total strain energy release rate (SERR) of multiscale models as a function of characteristic 

length. 
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